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Abstract: In this paper, wc consider a regression model built on depen- 
dent variables. This regression modelizes an input output relationship. Un- 
der boundedness type assumptions on the joint distribution function of the 
input variables, we show that a generalized Hoeffding-Sobol decomposi- 
tion is available. This leads to new indices measuring the sensitivity of the 
output with respect to the input variables. We also study and discuss the 
estimation of these new indices. 
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1. Introduction 

Sensitivity analysis (SA) aims to identify the variables that most contribute to 
the variability into a non linear regression model. Global SA is a stochatic ap- 
proach whose objective is to determine a global criterion based on the density 
of the joint probability distribution function of the output and the inputs of the 
regression model. The most usual quantification is the variance-based method, 
widely studied in SA literature. Hoeffding decomposition [9] (see also Owen 
[21]) states that the variance of the output can be uniquely decomposed into 
summands of increasing dimensions under orthogonality constraints. Following 
this approach, Sobol [26] introduces variability measures, the so called Sobol 
sensitivity indices. These indices quantify the contribution of each input on the 
system. 

Different methods have been exploited to estimate Sobol indices. The Monte 
Carlo algorithm was proposed by Sobol [27], and has been later improved by 
the Quasi Monte Carlo technique, performed by Owen [22]. FAST methods are 
also widely used to estimate Sobol indices. Introduced earlier by Cukier et al. [3] 
[4] , they are well known to reduce the computational cost of multidimensional 
integrals thanks to Fourier transformations. Later, Tarantola et al. [29] adapted 
the Random Balance Designs (RBD) to FAST method for SA (sec also recent 
advances on the subject by Tissot et al. [30]). 

However, these indices are constructed on the hypothesis that input variables 
are independent, which seems unrealistic for many real life phenomena. In the 
literature, only a few methods and estimation procedures have been proposed to 
handle models with dependent inputs. Several authors have proposed sampling 
techniques to compute marginal contribution of inputs to the outcome variance 
(see the introduction in Mara and references therein [17]). As underlined in Mara 
et aZ. [17], if inputs are not independent, the amount of the response variance 
due to a given factor may be influenced by its dependence to other inputs. 
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Therefore, classical Sobol indices and FAST approaches for dependent variables 
are difBeult to interpret (see, for example. Da Veiga's illustration [5] p. 133). Xu 
and Gertner [32] proposed to decompose the partial variance of an input into 
a correlated part and an uncorrelated one. Such an approach allows to exhibit 
inputs that have an impact on the output only through their strong correlation 
with other incomes. However, they only investigated linear models with linear 
dependences. 

Later, Li et al. [15] extended this approach to more general models, using the 
concept of High Dimensional Model Representation (HDMR [14]). HDMR is 
based on a hierarchy of component functions of increasing dimensions (trunca- 
tion of Sobol decomposition in the case of independent variables). The compo- 
nent functions are then approximated by expansions in terms of some suitable 
basis functions (e.g., polynomials, splines, ...). This meta-modeling approach 
allows the splitting of the response variance into a correlative contribution and 
a structural one of a set of inputs. Mara et al. [17] proposed to decorrelate the 
inputs with the Gram-Schmidt procedure, and then to perform the ANOVA- 
HDMR of Li et al. [15] on these new inputs. The obtained indices can be inter- 
preted as fully, partially correlated and independent contributions of the inputs 
to the output. Nevertheless, this method does not provide a unique orthogonal 
set of inputs as it depends on the order of the inputs in the original set. Thus, 
a large number of sets has to be generated for the interpretation of resulting 
indices. As a different approach, Borgonovo et al. [1. 2] initiated the construc- 
tion of a new generalized moment free sensitivity index. Based on geometrical 
consideration, these indices measure the shift area between the outcome density 
and this same density conditionally to a parameter. Thanks to the properties of 
these new indices, a methodology is given to obtain them analytically through 
test cases. 

Notice that none of these works has given an exact and unambiguous defi- 
nition of the functional ANOVA for correlated inputs as the one provided by 
Hoeffding-Sobol decomposition when inputs are independent. Consequently, the 
exact form of the model has neither been exploited to provide a general variance- 
based sensitivity measures in the dependent frame. 

In a pionnering work. Hooker [10], inspired by Stone [28], shed new lights on 
hierarchically orthogonal function decomposition. We revisit and extend here 
the work of Hooker. We obtain hierarchical functional decomposition under a 
general assumption on the inputs distribution. Furthermore, we also show the 
uniqueness of the decomposition leading to the definition of new sensitivity in- 
dices. Under suitable conditions on the joint distribution function of the input 
variables, we give a hierarchically orthogonal functional decomposition (HOFD) 
of the model. The summands of this decomposition are functions depending only 
on a subset of input variables and are hierarchically uncorrelated. This means 
that two of these components are orthogonal whenever all the variables involved 
in one of the summands also appear in the other. This decomposition leads to 
the construction of generalized sensitivity indices well tailored to perform global 
SA when the input variables are dependent. In the case of independent inputs, 
this decomposition is nothing more than the Hoeffding one. Furthermore, our 
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generalized sensitivity indices are in this case the classical Sobol ones. In the 
general case, the computation of the summands of the HOFD involves a mini- 
mization problem under constraints (see Proposition 1). A statistical procedure 
to approach the solution of this counstrainted optimization problem will be in- 
vestigated in a next paper. Here, we will focus on the particular case where the 
inputs are independent pairs of dependent variables (IPDV) . Firstly, in the sim- 
plest case of a single pair of dependent variables, the HOFD may be performed 
by solving a functional linear system of equations involving suitable projection 
operators (see Procedure 1). In the more general IPDV case, the HOFD is then 
obtained in two steps (see Procedure 2). The first step is a classical Hoeffding- 
Sobol decomposition of the output on the input pairs, as developped in Jacques 
et al. [11]. The second step is the HOFDs of all the pairs. In practical situ- 
ations, the non parametric regression function of the model is generally not 
exactly known. In this case, one can only have at hand some realizations of the 
model and have to estimate, with this information, the HOFD. Here, we study 
this statistical problem in the IPDV case. We build estimators of the general- 
ized sensitivity indices and study numerically their properties. One of the main 
conclusion is that the generalized indices have a total normalized sum. This is 
not true for classical Sobol indices in the frame of dependent variables. 
The paper is organized as follows. 

In Section 2, we give and discuss general results on the HOFD. The main result 
is Theorem 1. We show here that a HOFD is available under a boundedness type 
assumption (C.2) on the density of the joint distribution function of the inputs. 
Further, we introduce the generalized indices. In Section 3, we give examples 
of multivariate distributions to which Theorem 1 applies. We also state a suf- 
ficient condition for (C.2) and necessary and sufficient conditions in the IDPV 
case. Section 4 is devoted to the estimation procedures of the components of the 
HOFD and of the new sensitivity indices. Section 5 presents numerical applica- 
tions. Through three toy functions, we estimate generalized indices and compare 
their performances with the analytical values. In Section 6, we give conclusions 
and discuss future work. Technical proofs and further details are postponed to 
the Appendix. 

2. Generalized HoefFding decomposition- Application to SA 

To begin with, let introduce some notation. We briefly recall the usual functional 
ANOVA decomposition, and Sobol indices. We then state a generalization of this 
decomposition, allowing to deal with correlated inputs. 

2. 1 . Notation and first assumptions 

We denote by C the strict inclusion, that isAcB^AOBy^B, whereas we 
use C when equality is possible. 
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Let A, P) be a probability space and let Y be the output of a deterministic 
model rj. Suppose that 77 is a function of a random vector X = {Xi, • • • , Xp) € 
K^i P > 1 and that Px is the pushforward measure of P by X, 

{VL,A,P) {W,B{W),Px) (R,S(M)) 

Let 1/ he a. a-finite measure on (R^',S(Rp)). Assume that Px << v and let 

dPx 

"Px be the density of Px with respect to that is px ~ — ; — • 

av 

Also, assume that r\ G LKM^, ;B(K^), Px)- The associated innner product of this 
Hilbert space is: 

(/ll,/l2>= y /ll(x)/l2(x)pxrfi^(x) =E(/li(X)/l2(X)) 

Here E(-) denotes the expectation. The corresponding norm will be classically 
denoted by || • 1]. 

Further, F(-) = E[(- - E(-))^] denotes the variance, and Cov(-,*) = E[(- - 
E(-))(* — E(*))] the covariance. 

Let T'p := {1, • • • ,p} and S be the collection of all subsets of Vp. 
Define S~ :~ S\Vp as the collection of all subsets of Vp except Vp itself. 

Further, let X^ := (A';)ig„, u G S \ {0}. We introduce the subspaces of 
Ll(RP,B{W),Px) iHu)ueS, {H°)ueS and H° . iJ„ is the set of ah measurable 
and square integrable functions depending only on X„. TJg is the set of constants 
and is identical to {Hl^)ues- H^, w G 5 \ 0, and H° are defined as follows: 

H"^ = {KiXu) G i?„, (/i„, = 0, V « C M, V K e H°} 

H" = J h{X) = huiXu),K e \ 
[ ues J 

At this stage, we do not make assumptions on the support of X. For u G >S'\0, 
the support of X^ is denoted by ■ 



2.2. Sobol sensitivity indices 



In this section, we recall the classical Hoeffding-Sobol decomposition, and the 
Sobol sensitivity indices if the inputs are independent, that is when Px = 

Px^®■■■®Px^■ 

The usual presentation is done when X ~ W([0, 1]^) [26], but the Hoeffding de- 
composition remains true in general case [31]. 
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Let X = (xi, • • • , Xp) £ and assume that ij G L^(RP, Px)- The decomposi- 
tion consists in writting ?7(x) = r](xi, • • • , Xp) as the sum of increasing dimension 
functions: 



?7(x) = -no + ^r]i{xi) + ^ iji..j{xi,Xj) -\ |-?7i,...,p(x) 

i=l l<i<j<p 

VuiXu) (1) 

iiC{l---p} 

The expansion (1) exists and is unique under one of the hypothesis: 

i) J riu{xu)dPxi y i e u,y u C {!■ ■ -p} 
or 

ii) J r]uixu)r]vixy)dPx ^ V tt, v C {1 • • -p}, u^v 

Equation (1) tehs us that the model function Y ~ ?7(X) can be expanded 
in a functional ANOVA. The independence of the inputs and the orthogonality 
properties ensure the global variance decomposition of the output as V(Y) = 
EuesViVuiXu)). 

Moreover, by integration, each term has an explicit expression, given by: 



77o = E(X), 77, =E(y/X,)-E(y), ^ = ,p, 77,, = E(y/X„)-^ 77,, |u| > 2 

(2) 

Hence, the contribution of a group of variables X„ in the model can be 
quantified in the fluctuations of Y. The Sobol indices expressions are defined 
by: 

^""y(yy- v{Y) ' ""-^^ 

Furthermore, 



However, the main assumption is that the input parameters are independent. 
This is unrealistic in many cases. The use of expressions previously set up is 
not excluded in case of inputs' dependence, but they could lead to an unobvious 
and sometimes a wrong interpretation. Also, technics exploited to estimate them 
may mislead final results because most of them are built on the hypothesis of 
independence. For these reasons, the objective of the upcoming work is to show 
that the construction of sensitivity indices under dependence condition can be 
done into a mathematical frame. 
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In the next section, we propose a generalization of the Hoeffding decomposi- 
tion under suitable conditions on the joint distribution function of the inputs. 
This decomposition consists of summands of increasing dimension, like in Ho- 
effding one. But this time, the components are hierarchically orthogonal instead 
of being mutually orthogonal. The hierarchical orthogonality will be mathemat- 
ically defined further. Thus, the global variance of the output could be decom- 
posed as a sum of covariance terms depending on the summands of the HOFD. 
It leads to the construction of generalized sensitivity indices summed to 1 to 
perform well tailored SA in case of dependence. 

2.3. Generalized decomposition for dependent inputs 

We no more assume that Px is a product measure. Nevertheless, we assume: 



(C.l) 

= vildxi) (g) • • • (g) Vp{dxp) 
Our main assumption is : 

' (C.2) 




3 < M < 1, V u C Pp, px >M- px^Px^ 



where u'^ denotes the complement set of m in T^p. px^ and px^c are respec- 
tively the marginal densities of Xu and X„c . 

The section is organized as follows: a preliminary lemma gives the main result 
to show that is a complete space. Then, this ensures the existence and the 
uniqueness of the projection of ?/ onto . The generalized decomposition of rj 
is finally obtained by adding a residual term orthogonal to every summand, as 
suggested in [10]. The first part of the reasoning is mostly inspired by Stone's 
work [28], except that our assumptions are more general. Indeed, we have a 
relaxed condition on the inputs distribution function. Moreover, the support X 
of X is general. 

To begin with, let us state some definitions. In the usual ANOVA context, 
a model is said to be hierarchical if for every term involving some inputs, all 
lower-order terms involving a subset of these inputs also appear in the model. 
Correspondingly, a hierarchical collection T of subsets of "Pp is defined as follows: 

Definition 1. A collection T d S is hierarchical if for u € T and v a subset of 
u, one has w G T. 
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The next Lemma is a generalization of the Lemma 3.1 of [28]. As aheady 
mentioned, it will be the key to show the hierarchical decomposition. 

Lemma 1. Let T C S be hierarchical. Suppose that (C.l) and (C.2) hold. Set 
(5 = 1 — \/l — M g]0, 1]. Then, for any hu € H^, u G T, we have: 

E[(^ h^X)^] > nhli^)] (4) 

The proof of Lemma 1 is postponed to the Appendix. Our main theorem 
follows: 

Theorem 1. Let ?/ be any function in L'^{M.p , B{M.p), Px) . Then, under (C.l) 
and ( C.2), there exist functions rjQ^rji, - ■ ■ , rj-p^ G iJg x H'^ x • • • H!^ such that 
the following equality holds : 

Tj{Xi,--- ,Xp) = ^r7,(XO + 5]r;,,(X„X,) + ---+77P,(^i,--- 

i i,j 

= J2vu{Xu) (5) 

ues 

Moreover, this decomposition is unique. 

The proof is given in the Appendix. 
Notice that, in case where the input variables Xi, • • • , Xp are independent, S = I 
and Inequality (4) of Lemma 1 is an equality. Indeed, in this case, this equality 
is directly obtained by orthogonality of the summands. 

The variational counterpart of Theorem 1 is a minimization problem under 
conditional constraints. 

Proposition 1. Suppose that (C.l) and (C.2) hold. Let [V) he the minimization 
problem under constraints: 

r min - Y.^u{Xu)f] 

y E(7y„(X„)/X„\,) =0, V^eM, Vwg5\0 
Then (V) admits a unique solution rj* ~ {r]u)ues- 

Proof of Proposition 1 is postponed to the Appendix. Notice that a similar 
result for the Lebesgue measure is given in [10]. Its purpose was to provide 
diagnostics for high-dimensional functions. Here, we will no more exploit this 
idea. This will be done in a forthcoming work. Instead, we are going to construct 
stochastic sensitivity indices based on the new decomposition (5) and focus on 
a specific estimation method for IPDV models. 
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2.4- Generalized sensitivity indices 

As stated in Theorem 1, under (C.l) and (C.2), the output Y of the model can 
be uniquely decomposed as a sum of hierarchically orthogonal terms. Thus, the 
global variance has a simplified decomposition into a sum of covariance terms. 
So, we can define generalized sensitivity indices. 

Definition 2. The sensitivity index Su of order \u\ measuring the contribution 
of Xu into the model is given by : 



y(77„(x„)) + E„ 



Cov{riuiXu),r]y{X.„)) 



V{Y) 



(6) 



More specifically, the first order sensitivity index St is given by : 



V{v^{X,))+J:u^<iCov{7^,{X,),r,,{X^)) 



An immediate consequence is given in Proposition 2 (sec proof in the Ap- 
pendix) : 

Proposition 2. Under (C.l) and (C.2), the sensitivity indices Su previously 
defined are summed to 1, i.e. 

«es\{0} 

Thus, sensitivity indices arc summed to 1. Furthermore, the covariance terms 
included in these new indices allow to take into account the inputs dependence. 
Thus, we arc now able to measure the influence of a variable on the model, espe- 
cially when a part of its variability is embedded into the one of other dependent 
terms. We can distinguish the full contribution of a variable and its contribution 
into another correlated income. 



(7) 



Note that for independent inputs, the summands rju are mutually orthogonal, 
so Cov(77„, rju) = 0, u ^ V, and we recover the well known Sobol indices. Hence, 
these new sensitivity indices can be seen as a generalization of Sobol indices. 

However, the HOFD and subsequent indices arc only obtained under con- 
straints (C.l) and (C.2). In the following, we give illustrations of distribution 
functions satisfying these main assumptions. 
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3. Examples of distribution function 

This section is devoted to examples of distribution function satisfying (C.l) and 
(C.2). The first hypothesis only implies that the reference measure is a product 
of measures, whereas the second is trickier to obtain. 

In the first part, we give a sufficient condition to get (C.2) for any number p 
of input variables. The second part deals with the case p = 2, for which wc give 
equivalences of (C.2) in terms of copulas. 

3.1. Boundedness of the inputs density function 

The difficulty of Condition (C.2) is that the inequality has to be true for any 
splitting of the set {Xi, ■ ■ ■ ,Xp) into two disjoint blocks. We give a sufficient 
condition for (C.2) to hold in Proposition 3 (the proof is postponed to the 
Appendix) : 

Proposition 3. Assume that there exist Mi,M2 > with 



Ml <Px< M2 (C.3) 



Then, Condition (C.2) holds. 

Let give now an example where (C.3) is satisfied. 

Example 1: Let v be the multidimensional gaussian distribution Np(m, E) 
with 

rl ... 





/ mA 




m = 








\mp) 





... 

Assume that Px is a Gaussian mixture a ■ Np{ni,J^) + {1 — a) ■ Np{fi,n), 
a elO,l[ with 



/vl Pi2 ■■■ Pip' 




kP'p/ 

Then, (C.3) holds iff the matrix (57^^ — S^^) is positive definite. 

In the next section, we will see that (C.2) has a copula version when p — 2. 
We will give some examples of distribution satisfying one of these conditions. 
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3.2. Examples of distribution of two inputs 

Here, we consider the simpler case of inputs X = {Xi,X2). Also, until Section 
4, we will assume that v is absolutely continuous with respect to Lebesgue mea- 
sure. The structure of dependence of Xi and X2 can be modelized by copulas. 
Copulas [19] give a relationship between a joint distribution and its marginals. 
Sklar's theorem [25] ensures that for any distribution function F{xi^X2) with 
marginal distributions _F'i(a;i) and F2{x2), F has the copula representation, 

F{xi,X2) = CiF,{x,),F2ix2)) 

where the measurable function C is unique whenever Fi and F2 are absolutely 
continuous. 

The next corollary gives in the absolutely continuous case the relationship 
between a joint density and its marginal: 

Corollary 1. In terms of copulas, the joint density ofX. is given by: 

Px{xi,X2) = c{Fi{xi),F2{x2))pxi{xi)pxAx2) (9) 



Furthermore, 



c{u, v) 



dudv 



(10) 



Now, Condition (C.2) may be rephrased in terms of copulas: 



Proposition 4. For a two-dimensional model, the three following conditions 
are equivalent: 



px > M ■ PX1PX2 v-a.e. for some < M < 1 (C.4) 



2. c{u, v)>M, V {u, v) e [0, 1] 



(C.5) 



C{u,v) — Muv + (1 — M)C{u,v), for some copula C (C.6) 



The proof of Proposition 4 is postponed to the Appendix. 
Hence, the generalized Hoeffding decomposition holds for a wide class of 
examples. The first example is the Morgenstern copulas [18]: 

Example 2: The expression of the Morgenstern copulas is given by: 

Coin, v) = uv[l + 61(1 - u){l - v)], 9 e [-1, 1] 
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For 6* e] - 1, 1[, (C.6) holds, and 

Ceiu, v) = {l- \9\)uv + \9\uv[l + u){l - v)] 

Let now consider the class of Archimedian copulas, 

C{u,v) = ip-'[^{u) + ^{v)], u,«e[0,l] (11) 

where the generator ip is sl non negative two times differentiable function defined 
on [0, 1] with ip{l) = 0, tp'iu) < and ip"{u) > 0, V u S [0, 1]. 

A sufficient condition for (C.5) is given in Proposition 5: 
Proposition 5. // there exist AIi, M2 > such that: 

- ip'{u) > A/i V u e [0, 1] (12) 
^il^^)>M,,yue[0,l] (13) 

Then, Condition (C.5) holds. 

The proof is straightforward. Now, we will see three illustrative Archimedian 
copulas satisfying (C.5). 

Example 3: The Frank copula is characterized by the generator: 

^i(x)=log(^^^-^), (?eR\{0} 

Condition (C.5) holds, and c(u,v) > -e{er^ - l)e-2« if 6* > 0, c(w, w) > 
—9(e~^ — 1) elsewhere. 

The next two examples also satisfy (C.5) by the intermediate Proposition 5. 
Example 4-' Let a < 0, 9 > and /3 with f3 < —ae~^. Set 

^2(^)--?e-'^^+/3x + (^e-''-/3), xe[0,l] (14) 



Example 5: Let C < and set 

ip:i{x)=x\nx + {C-l)x + {l~C), x S [0, 1] (15) 



Leaving the class of copulas, we now directly work with the joint density 
function. Proposition 6 gives a general form of distribution for our framework: 
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Proposition 6. If px has the form 

Pxixi,X2) ^ a ■ fxA^i)fx2{x2) + - a) ■ 9x{xi,X2), Q!G]0, 1[ (16) 

where fx^ , fx2 o,'"'^ univariate density functions, and gx is any density function 
(with respect to v) with marginals fxi and fx^i then px satisfies (C.5). 

The proof is straightforward. 

Example 6: As an ihustration of Proposition 6, take v = vl, fx = fxifx2 
a normal density with a diagonal covariance matrix E, and gx a normal density 
of covariance matrix fl, with flu = E^^, i = 1,2. Notice that because a copula of 
Gaussian mixture distribution is a mixture of Gaussian copulas (see [20]), this 
example can be directly recovered by the copula approach. 

Example 7: Let generalize Example 6. If Px is a Gaussian mixture 
a- N2{m,Y.) + {l-a)- N2{tJ.,n), ae]0,l[ 

with 

then (C.4) holds iff uj^ < af and w| < fT|. 

Thus, for many distributions, the generalized decomposition holds, and gen- 
eralized sensitivity indices may thus be defined. 

For the remaining part of the paper, we will assume that the set of inputs is 
an IPDV. If p is odd, we will assume that an input variable is independent to 
all the others. 

The next section is devoted to the estimation of HOFD components. The sim- 
plest case of a single pair of dependent variables is first discussed. Then, the 
more general IPDV case is studied. In this last part, first and second order 
indices are defined to measure the contribution of each pair of dependent vari- 
ables and each of its components in the model. Indices of order greater than one 
involving variables from different pairs will not not be studied here. 

4. Estimation 

Using the property of hierarchical orthogonality (iJ^ _L H^, V v C u), we will 
see that the summands of the decomposition arc solution of a functional linear 
system. For u G S, the projection operator onto iJ" is denoted by Pj^o . 

In this section, we present the HOFD terms computation, based on the res- 
olution of a functional linear system. The result relies on projection operators 
previously set up. Further, we expose the linear system estimation for practice. 
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4-1. Models of p — 2 input variables 

This part is devoted to the simple case of bidimensional models. Let 

r = 77(Xi,X2) 

Assuming that Conditions (C.l) and (C.2) both hold, we proceed as follows: 
Procedure 1 



1. HOFD of the output: 

r = 7?o + in{Xi) + i^2{X2) + m2{Xi,X2) (17) 

2. Projection of Y = t]{X) on H°, V it C {1, 2}. As H° ± H^,\f v C u, we 



obtain: 



fid 
Id 









Pho 

Id 

IdJ 



fvo\ 



m 
\1n2J 



Phi 
\0 

3. Computation of the right-hand side vector of (18): 

( 



' Ph«{v) * 

\PhiMJ 



(18) 



(PhM\ 

Pho{v) 
PHoiv) 

\PhoJv)J 



nv/Xi)-Eirj) 
V/ - E(?7/^i) - E(7?/X2) + E(ry) / 



(19) 



In this frame, we have: 

Proposition 7. Letr] be any function of L'^{Rp , B{W) , Px) . Then, under 
(C.l) and (C.2), the linear system 



TrI I 



(Id 



Vo 



Ph« 

Id 









admits in h = (/iq, 





Id 



IdJ 

hi2) e Hid X 



/ho\ 
hi 
h2 
\hr2j 



(PhM\ 

PH«iv) 
PHoiv) 
\PhoJv)J 

X H12 the unique solution h 



(20) 



Reduction of the system (18). As the constant term corresponds to the 
expected value of rj, and the residual one can be deduced from the others, 
the dimension of the system (20) can even be reduced to: 



Id 
Ph'> 



Pro 
Id 



E{j]/Xi)-E{rj) 
E(?7/X2)-E(r/) 



(21) 
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5. Practical resolution: The numerical resolution of (21) is achieved by an 
iterative Gauss Seidel algorithm [13] which consists first in decomposing 
A2 as a sum of lower triangular (L2) and strictly upper triangular (U2) 
matrices. 

Further, the technique uses an iterative scheme to compute A. At step 
fc + 1, we have: 



A(^+^) := (^[.+1)1 =i2-^(B-C/2- AW) (22) 
Using expression of A2, we get: 

ny- A^'V^i) - nr a^) \ 



6. Convergence of the algorithm: now, we hope that the Gauss Seidel algo- 
rithm converges to the true solution. Looking back at (18) , we see that we 
only have to consider Pfja ( respectively Pfjo ) restricted to H2 ( respectively 
toH^). 

Under this restriction, let us define the associated norm operator as : 
\\PHof sup E[Pho{U)% ^,J = 1,2, j ^ i 

E{U^) = l 

As explained in [7], Gauss Seidel algorithm converges to the true solution 
A if A2 is striclty diagonally dominant, which is implied by : 

||Pffo|| < t = l,2 (24) 

i.e. 

sup E[P^o(C/)2] < 1, i^l,2, j ^i (25) 

E(C/^) = l 

Equality (19) reveals that Ph«{U) = E{U/X,) - E{U). Hence, by the 
Jensen inequality [12] for conditional expectations, \\P}jo\\,i = 1,2 ad- 
mits an upper bound: 



Take U e 



\Pho\\ = sup E[(E(C//X2) - E([/))2] 

E(!7^) = l 

ueH^ 

sup E[E({//X2)^] asU e H[ 
< sup E[E(C/V^2)] = 1 by Jensen 

E((7^) = l 
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The same holds for U G H2, and we also have ||^^fp|| < 1. 
Moreover, the Jensen's inequality is strict ifU is not Xi-measurable. As U 
is a function of Xj (that is j = 2 if i = 1 and conversely), the condition of 
convergence holds if Xi is not a measurable function of X2 ■ Hence, Gauss 
Seidel algorithm converges if Xi is not a function of X2 ■ 

7. Estimation procedure: Suppose that we get a sample of n observations 

• Estimation of the components of the HOED: the iterative scheme 
(23) requires to estimate conditional expectations. As extended in Da 
Veiga et al [6], we propose to estimate them by local polynomial re- 
gression at each point of observation. Then, we use the leave-one-out 
technique to set the learning sample and the test sample. Moreover, 
as the local polynomial method can be summed up to a generalized 
least squares (see Ean and Gijbels [8]), the Sherman- Morrison for- 
mula [24] is applied to reduce the computational time. 

A more detailed procedure is given in the Appendix. The iterative al- 
gorithm is easy to implement. We stop when IjA'-'^^"'^' — A'*"'^!! < e, 
for a small positive e. 

Once (771,772) have been estimated, we estimate rjo by the empirical 
mean of the output. Then, an estimation of r;i2 is obtained by sub- 
straction. 

• We use empirical variance and covariance estimation to estimate sen- 
sitivity indices Si, S2 and S12. 

4-2. Generalized IPDV models 

Assume that the number of mputs is even, so p = 2k, k > 2. We note each group 
of dependent variables as X*^'^ := {x[^\ X2^), i ~ 1, ■ ■ ■ ,k. By rearrangement, 
we may assume that: 

X = {Xi,X2, ■ ■ ■ , X2k-l,X2k) 

SA for IPDV models has already been treated in [11]. Indeed, they proposed 
to estimate usual sensitivity indices on groups of variables via Monte Carlo 
estimation. Thus, they have interpreted the influence of every group of variables 
on the global variance. Here, we will go further by trying to measure the influence 
of each variable on the output, but also the effets of the independent pairs. 
To begin with, as a slight generalization of [26] and used in [11], let apply the 
Sobol decomposition on groups of dependent variables, 

fc 

\u\=2 
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where for u = {ui,--- ,ut} and t = we set X^") = {X'^^^\--- 
Furthermore, (?7„, r/^) = 0, V u 7^ w. 

Thus, we obtain independent components of IPDV. Under the assumptions 
discussed in the previous section, we can apply the HOFD of each of these 
components, that is, 

with {ipi^u, ^i,v) = 0, V w C u C {1,2}. In this way, let define some new 
generalized indices for IPDV models: 

Definition 3. For i ~ 1, ■ ■ • ,k, the first order sensitivity index measuring the 
contribution of x[^'^ (respectively ) on the output of the model is: 



r, V{ip,A + Cov{Lpi,l,Lpi^2) _V{ip,^2) + Cov{(p,,l,Vt.2)\ 

5m- ' (^5^2- j (26) 

The second order sensitivity index for the pair X^^\ i — 1, - ■ ■ ,k, is defined 

as: 

602 - (27) 

The estimation procedure of these indices is quite similar to Procedure 1: 
Procedure 2 

1. Estimation of (''7i)i=i,- - ,fe' as reminded in Part 2.2 with Equations (2), 
7}i = E(l'/X(*)) — E(y). Step 7 of Procedure 1 gives method to estimate 
the conditional expectations. So that, we will have estimations of r]i, i = 
I,--- ,fc. 

2. For i ~ I, ■ ■ ■ , k, we apply step 2 to step 7 of Procedure 1, considering rji 
as the output. 

If p is odd, the procedure is the same except that the influence of the inde- 
pendent variable is measured by a Sobol index, as it is independent from all the 
others. 

The next part is devoted to numerical examples. 



5. Numerical examples 

In this section, we study three examples with dependent input variables. We 
consider IPDV models and a Gaussian mixture distribution on the input vari- 
ables. We choose covariance matrices of the mixture satisfying conditions of 
Example 1. 
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We give estimations of our new indices, and compare them to the analytical 
ones, computed from expressions (6). We also compute dispersions of the esti- 
mated new indices. 

In [6], Da Veiga et al. proposed to estimate the classical Sobol indices Su = 
{V[E{Y/Xu)] - J2v(iu V[E(Y/Xy)])/V{Y),uC Vp, by nonparametric tools. In- 
deed, the local polynomial regression were used to estimate conditional moments 
E{Y/Xu), u C Vp. This method, used further, will be called Da Veiga procedure 
(DVP). Results given by DVP arc compared with the ones given by our method. 
The goal is to show that the usual sensitivity indices are not appropriate in the 
dependence frame, even if a relevant estimation method is used. 

5.1. Two-dimensional IP DV model 

Let consider the model 

Y ^Xi+X2 + X1X2 
Here, v and Px are of the form given by Example 1. with m = = 0. 

Thus, the analytical decomposition of Y is 

7?o = E(XiX2), ryi = Xi, 7]2 = X2 r?i2 = X1X2 - E{XiX2) 

For the application, we implement Procedure 1 in Matlab software. We pro- 
ceed to L = 50 simulations and n ~ 1000 observations. Parameters were fixed 
at (Ti = (72 = 1, 'Pi=f2= 0-5> P12 ~ 0.4 and a = 0.2. 

In Table 1, we give the estimation of our indices and their standard deviation 
(indicated by ±-) on L simulations. In comparison, we give the analytical value 
of each index. 

The analytical classical Sobol indices are difficult to obtain, but we give 
estimators of the classical Sobol indices with DVP. 



Table 1 

Estimation of the new and DVP indices with pi2 = 0.4 







5i 


S2 


5l2 




New 


Estimation 


0.42 ± 0.041 


0.41 ±0.043 


0.17 ±0.026 


l±9.10-is 


indices 


Analytical 


0.39 


0.39 


0.22 


1 


DVP 
indices 


Estimation 


0.64 ±0.045 


0.65 ±0.044 


0.41 ± 0.038 


1.7 ±0.09 



We notice that estimations with our method give quite good results in com- 
parison with their analytical values. The estimation error of the interaction term 
is due to the fact that the component 7)12 is obtained by difference between the 
output and the other estimated components. 
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The DVP indices are arc difficult to interpret as the sum is higher than 1. 

In our method, it would be relevant to separate the variance part to the 
covariance one in the first order indices. Indeed, in this way, we would be able 
to get the part of variability explained by Xi alone in Si, and its contribution 
hidden in the dependence with Xj. We note the variance contribution alone, 
and Sf the covariance contribution, that is 



Coy{X„Xj 




1,2, j^i 



The new indices estimations given in Table 1 are decomposed in Table 2. As 
previously, the number at the right of ± indicates the standard deviation on L 
simulations. 



Table 2 

Estimation of 5" and S'r with pi2 = 0.4 





St 




Si 




0.28 ±0.04 


0.14 ±0.01 


0.42 ± 0.041 




0.27 ± 0.043 


0.14 ±0.01 


0.41 ± 0.043 


Analytical 


0.25 


0.14 


0.39 



For each index, the covariate itself explains 28% (in estimation, 25% in re- 
ality) of the part of the total variability. However, the contribution embedded 
in the correlation is not negligible as it represents 14% of the total variance. 
Considering the shape of the model, and coefficients of parameters distribution, 
it is quite natural to get the same contribution of Xi and X2 into the global 
variance. Also, as their dependence is quite important with a covariance term 
equals to 0.4, we are not surprised by the relatively high value of SI (resp. S"!). 

From now, we take pi2 = 0, i.e. we assume that the inputs are independent. 
Let compare our new estimated indices with their analytical values in Table 3. 
We again decompose new indices into a variance {S^) and a covariance (Sf) 
contribution. 



Table 3 

Comparison between analytical and estimated indices with p\2 = 





New indices 


Theoretical 






s. 


s. 




0.39 ±0.039 


-0.01 ±0.01 


0.38 ±0.036 


0.375 


X2 


0.38 ± 0.045 


-0.01 ±0.024 


0.37 ±0.037 


0.375 


XIX2 


0.26 ±0.038 


-0.01 ±0.01 


0.25 ±0.027 


0.25 



Thus, the new indices are well tailored if we have a small idea on inputs depen- 
dence in a system. Indeed, Table 3 shows that our new indices take dependence 
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into account if it exists, and the covariance contribution is estimated by if not. 
New indices recover the classical Sobol indices in case of independence. 

5.2. Linear Four- dimensional model 

The test model is 

Y = 5Xi+ 4X2 + 3X3 + 2X4 

Actually, Condition (C.2) only needs to be satisfied on groups of correlated 
variables. 

Let consider the two blocks X*^^^ = (XijX^) and X'-^^ = {X2,X4^) of correlated 
variables. 

The previous form of density can be taken for X'^' and X'-^K Px'-') then the 
Gaussian mixture a,; • iV2(0,l2) + (1 — a^) • A^2(0,i^i), i = 1,2. The analytical 
sensitivity indices are given by (26) & (27). 

2(1) 2(1) 

For L = 50 simulations and n = 1000 observations, we took ip^ ~ = 
0.5, (^1^^^ = 0.7, v?2^^^ = 0.3, pil = 0.4, pyl = 0.37 and ai = aa = 0.2. 

Figure 1 displays the dispersion of indices of first order for all variables and 
second order for grouped variables. We compare them to their analytical values. 
In the same figure, we also represented the estimators of classical Sobol indices 
with DVP. 
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* Analytical Index 
■ DVP Index 




o 



i 

S3 S4 SI 3 

Sensitivity Indices 



Fig 1. Boxplots representation of new indices- Comparison with analytical and DVP indices 



We see that Xi has the biggest contribution, whereas the influence of X4 
is very low. It reflects well the model if we look at the coefficients Xi, i = 
1, • • • ,4. Also, interaction terms are well estimated, as they are closed to 0. For 
each case, the dispersion on 50 simulations is very low. 

As for the DVP estimation, it is once again very high compared with the true 
indices values. 



5.3. The Ishigami function 

This function is well known in SA ([23]). It is defined by: 

Y ^ sm{Xi) + asm^{X2) + bX^ sm{Xi) 

We assume that ^"3 is the independent variable, and that Xi and X2 are corre- 
lated. Px is again the Gaussian mixture a ■ N3{0, 13) + (1 — a) • A^3(0, fi). 

With L = 50 simulations of n = 1000 observations, we fixed parameters of 
distribution at ipj = 0.15, ipl = 0.85, (fil = 0.75, pi2 = 0.3 and a = 0.2. 
In Figure 2, the dispersion of the new measures is represented for fixed 6 = 0.1 
and different values of a. In addition, the analytical new indices and estimated 
classical Sobol indices with DVP are displayed. 
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a= 0.1 a= 2.1 




SI S2 S3 S12 SI S2 S3 S12 

Sensitivity indices Sensitivity indices 



Fig 2. Boxplots of new indices for different values of a for Ishigami function 

Every boxplot shows that there is a small dispersion. For all estimations, 
DVP indices are larger that the new ones. The figure clearly shows that, for all 
values of a, the sum of these four indices is greater than 1. It again shows their 
non adaptation to a situation of dependence. 

If we have a look on the values taken by our new sensitivity indices, we see 
that, for small values of a, the variable Xi contributes the most to the model's 
variability. This role decreases as a increases, and X2 then gets the biggest 
contribution. For any value of a, the input plays a very negligible role, which 
seems realistic as 6 is a small fixed value. As for the interaction index S'12, it is 
getting bigger with the increasing importance of a, but the contribution remains 
low. 

6. Conclusions and Perspectives 

The Hocffding decomposition and associated Sobol indices have been widely 
studied in SA over past years. Recently, a literature appears to treat the case of 
dependent input variables, in which authors propose different ways to deal with 
dependence. The goal of this paper is to conciliate the problem of inputs depen- 
dence with the Hoeffding decomposition. Indeed, we study a functional ANOVA 
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decomposition in a generalized inputs distribution frame. Thus, we show that a 
model can be uniquely decomposed as a sum of hierarchically orthogonal func- 
tions of increasing dimension. Also, this approach generalizes the HoefFding's 
one, as we recover it in case of independence. 

Similarly to the classical Sobol decomposition, this leads to the construction of 
new sensitivity indices. They consist of a variance and a covariance contribution 
able to take into account the possible correlation among variables. In case of 
independence, these indices are the classical Sobol indices. However, the indices 
construction is only possible under specific assumptions on the joint distribution 
function of the inputs. We expose few cases that satisfy these assumptions for 
any p-dimensional models. More specifically, for two-dimensional models, the 
required assumption is equivalent to assumptions on copulas. In this context, 
we give examples satisfying one of these assumptions . 

Focused on the IPDV models, summands of the decomposition are estimated 
thanks to projection operations. This leads to the numerical resolution of func- 
tional linear systems. The strength of this method is that it does not require to 
make assumptions on the form of the model or on the structure of dependence. 
We neither use meta-modelling and avoid in this way many sources of errors. 
Through three applications on test models, we observe the importance of con- 
sidering the inputs correlation, and show how our method could catch it. The 
comparison with estimators of classical indices with DVP shows that the Sobol 
indices are not appropriate in case of correlations, even when using nonpara- 
metric method. Also, when inputs independence holds, the new indices remain 
well suited to measure sensitivity into a model. 

Nevertheless, only considering IPDV models for estimation is restrictive. The 
perspective is to explore other estimation methods suitable for more general 
models. Also, we intend to lead a systematic study on copulas satisfying or not 
our assumptions. 

Appendix A: Generalized Hoeffding decomposition 
A.l. Generalized decomposition for dependent inputs 

The upcoming proof follows the guideline of the proof of Lemma 3. 1 in Stone [28]. 
Proof of Lemma 1 

By induction on the cardinal of T, let show that 

H{n) : V r/#(T) = n, E[(E„eT ^"(X))^] > E„eT ^[^^(X)] 

• IS obviously true, as T is reduced to a singleton 

• Let 71 £ N*. Suppose that T-iln) is true for all 1 < n < n. Let T such that 
#{T) — n + 1. We want to prove 'H{n + 1). 

Choose a maximal set r ofT, i.e. r is not a proper subset of any set u inT. We 
show first that 
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E[(^/i„(X))'] > M-E(/i?(X)) (28) 

- U#{r) =P, by definition of H° , we get E[(X;„gT /i«(X))2] > E{hl{X)) > 
ME{hl{X)) as M <1. 

-If^< < P - 1, set X = (Xi.Xa), where Xr = {Xi)i^r and X2 = 

(X;);gr- By Condition (C.2), it follows that 

px > M ■ pxiPX2 

As a consequence, 



mJlueThuC^)?] = J J_^^[hr{x2) + Y:u^^h4xi,X2)]''pxHdxi,dx2) 



> M Jxjx2^^^^^'^'> + j2u^r^^(^^^^2)]'^PXiPX2^l{dXi)v2{dX2) 

> M J^^E[{hr{X2)+J2ujtr'^uixi,X2)f]px,l^lidXl) 



By maximality of r and by definition of H^, 

* If u C r , hu only depends on X2 and by orthogonality, 

E{h4X2)hr{X2)) ^ 

* If u <^ r, hu depends on Xi fixed at xi, and X2 ~ {Xi)i^rnu, so 
hu £ Hrnu! with r Hu C r, it comes then 

E{hu{xi,X2)hr{X2))^0 

Thus, 

E[(^/i„(X))2] > uf E{hl{X2))px,i^iidxi) 
= M-E{hriX)) 

So (28) holds for any size of any maximal sets of T. 

By using (28) with hr ~ hr and hu = — \l u ^ r, we get 



E[(/i,(X) -PYI M^)f] > ME{hl{X)) (29) 

E[/i,.(X)5] fe„(X)] 
Taking P = " ^^-.^oi > it follows that: 

E[(fe,(X)-/3E„^./*4X))'] > ME{hUy^)) 

that IS E[K{X)] huiX))^] - ^'^^C^ri^)) 

Hence, 
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[E{hriX) J2 huiy.))f < (1 - M) ■ E(/i?(X)) • E[(^ h4X)f] (30) 



This implies 



E[(^ hu{x)f] > (1 - VT^) 

u 

Set X — fer(X) and y = hs(X.).(31) is rephrased as 



\\x + yf>{l-Vl^){\\x\\^ + \\y\\'} 
Further, (30) is {x,y} > -VI - A/||2;|| ■ ||y||. Thus, 



\\xf + \\yf > 2{x,y} 
> 



(a^,j/) by (32) 



(31) 



(32) 



H(7i) is supposed to be true and (31) holds, it follows that: 

niT.uh.i^)?] > <5[E(/i?(X)) + <5"-^E.^.E(ft^(X))] 
> '5"E„E(/iS(X)) as Se ]0,1] 
= 5#(^)-iE„E(h^(X)) 

Hence, 'H{n + 1) holds. 
We can deduce that T-L{n) is true for any collection T ofVp. m 
Proof of Theorem 1 

Let define the vector space 74"" — {Eugs- hu{Xu), hu € H^u-^ G 'S' }• 

In the first step, we will prove that K° is a complete space to prove the existence 
and uniqueness of the projection of rj in , thanks to the projection theorem [16]. 
Secondly, we will show that rj is exactly equal to the decomposition into , and finally, 
we will see that each term of the summand is unique. 

• We show that H° is closed into Hu (as Hu is a Hilbert space). 

Let {h„,u)n be a convergent sequence of H'^ with h„,u hu. As {h„^u)n £ C 
Hu complete, hu € Hu . Let v C u, and hv € H° : 



{hu — h hu) = {hu, hu) —{h hu) 




as m ± H, 
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Thus, {hu, hy) — 0, so that hu G H^. is then a complete space. 

Let {hn)„ be a Cauchy sequence in K° and we show that each component is of 
Cauchy and that h„ ^ h £ . 

As h„ £ K'^ , hn = X^ugs- hn,u, hn.u G H^u- H follows that : 

> §*{S \\h„,u - hm,u\\'^ hy Inequality (4) 

As {hn)n is a Cauchy sequence, by the above inequality, {hn,u)n is also Cauchy. 
As hn u ^ hu € H^, we deduce that hn — > ^a- hu = h £ i^". 

Thus, K" is complete. By the projection theorem, we can deduce there exists a 
unique element into such that : 

• Decomposition ofrj: following Hooker [10], we introduce the residual term as 

tigS- 

By projection, (rj — Ylves- = V u G S~ , \/ hu € Hu- Hence, ri{X.) = 

Sugs VuiXu), rju G Hu, V It G S, and this decomposition is well defined. 

• Terms of the summand are unique: assume that rj = X^ues ~ X^ues ^ 
Hu- 

By Lemma 1, it follows that 



Proof of Proposition 1 

Let first prove the following equivalence 



J r]u{xu)r]v{xv)px{x)dv{x) = Vu C m, V r]v 
t 

J Tlu{xu)px{x)dvi{xi) duu^ixu") = Vu G S, Vi G M 



Let V C u and i £ u\v, then 



J r]u{xu)riu{xu)px{x)du(x) = / r)u{xu)riu{xu)px{,x)dvi{xi) dvu^ixu") di'u\i{xu\i) 

r)v(Xv) {^j ■nu{Xu)px{x)dVr{Xi) dUu^iXu")^ dVu\i[Xu\i) 



by assumption 



imsart-ejs ver. 2011/12/06 file: ps-template.tex date: March 13, 2012 



G. Chastaing et al. /Generalized ANOVA Decomposition for Dependent Variables 26 

Conversely, assume that 3 u, 3i £ u such that J rju{xu)px{x)dvi{xi) dj^„c(x„c) 7^ 
0, then, 

rjv — J riu{xn)px{x)dvi{xi) dvu^ixu") with v = u\i 



and 



j r]u{xu)rjv{xv)px{x)di/{x) = J r]u{xu) (^J r)u{xu)px{x)dvi dvu<^ px{x)dvi dv^c dv^\i 

riu{xy.)px{x)du^{x^) di/uc(a;„c)^ dv^\^{x^\i) 



> 

There is a contradiction, so that J r]u{xu)px{x)dvi{xi) dvu<:{xu<--) = Vi G u,\/u. 
The second expression can he rewritten as : 

j r)u{xu)px{x)dvi{xi) dvu<=ixu'=) = E(77„/X„\i) Vi G u.Vu G S 

Then, by Theorem 1, the minimization problem (V) admits a unique solution. 



A. 2. Generalized sensitivity indices 

Proof of Proposition 2 

Under (C.l) and (C.2), Theorem 1 states the existence and the uniqueness decom- 
position of Tj: 

77(X) = ^r,4X„), 

ties 

with H° ± H°, \f V Cu. Therefore, 



and 



V{Y) = VMX)) = E{r,\X)) - r,^o 

= J2^il^Ux^)) + J2nnu{XMx.)) 

= Y.'^{nu{x^))+Y. E nnu{x^),n.{x^)) 



ti/0 



■117^0 U7^0 
u^Vp ti^t),t)0ti 



E 

117^0 



V(j7t,(X„))+ Cov{r,n{X^),r,4X^)) 

V 

Thus, (6) holds, and equalities (7) and (8) follow obviously. 
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Appendix B: Examples of distribution function 
B.l. Boundedness of the inputs density function 

Proof of Proposition 3 

Let u C Vp, and < Mi < px < M2- As px„ and px^ are marginals, they are upper 
bounded by M2 . 

m| Mi _2 

As a consequence, px^px^c < Tir^i — HtPx, so that px > M1M2 Px^Px^c , 
^ All All " 

with < AfiAff^ < 1. 

■ 

Proof of Example 1 

• v is a product of measure as — — = 11^=1 ^ii^i)} with Vi ~ N(mi, aj) . So 



V 




1=1 



• Px is given by 



dPx , . dPx dvL , , 



1/2 



exp — i*(x — 7Ti)(Sl ^ 



S"')(x-m) (33) 



First, we have px(x) > q > 0. 

Further, the sufficient and necessary condition to have px < Ah < <x is to get 
{Q^^ — S^^) positive definite. Indeed, if {Q,^^ — S^^) admits a negative eigen- 
value, Px can not be bounded. Thus, < a < px < A/2 iff (O^^ — E~^) is 
positive definite. 



B.2. Examples of distribution of two inputs 

Proof of Proposition 4 

Condition (C.5) is immediate with Equation 9. Let prove that (C.5) is equivalent 
to (C.6). 

If (C.6) holds, then c(u,v) > AI . Conversely, we assume that < M <1, and 

=, , , _ C{u, v) ~ AIuv 
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It is enough to show it is a copula : Obviously, C{0,u) — C{u,0) = and C{l,u) — 

" ciu V) — I 

C{u,l) = u V It G [0, 1] . By second order derivation, it comes that c(u, v) = — — — ^ 

so c{u,v) > by hypothesis (C.5). 



Appendix C: Estimation 

C. 1 . Model of p = 2 input variables 

Proof of Proposition 7 



• We first show first that (S ) admits an unique solution. 

Under (C.l) and (C.2), by Theorem 1, the decomposition of 77(X) is unique and 
ri{Xi,X2) = Vo + Vi{Xi) + riiiXi) + Vi2{Xi, X2) 

with 

Vo e Hid 

rj,eH°±H@, i = l,2 
7712 G -L Hf, i = l,2, H°2 -L 



Thus, 



/Id 








'\ 





Id 










PjjO 

"2 


Id 





\o 








Id/ 



So {rjo , rji , ri2 , r]i2) is solution of (S). 



m 

\V12/ 



Pho{v) 
Pnoiv) 
\PhoM/ 



(34) 



Now, assume there exists an another solution of the system, say {r]o,- ■ ■ ,iTPp) G 
Hijj X ■ ■ ■ X , then 



'70 - J?o = 

rii —rji + Pjjo {ri2 — rj2) =0 
Pho (j?! - 571) + '72 - f?2 = 
, ?7i2 — 5712 = 



r]o = ^70 

Pn^lim - m + 172 - '72) 
Pn^iVi - f?i + 172 - '72) 
I, ?7i2 = rji2 



m - 571^+ V2-m& n Hi 

r]i2 = '012 

As i]i — 771 G Hi and r;2 — 772 G H2, it follows that 



imsart-ejs ver. 2011/12/06 file: ps-template.tex date: March 13, 2012 



G. Chastaing et al. /Generalized ANOVA Decomposition for Dependent Variables 29 



r (t?! — ^1 — J]! + ??2 — 772) = ^ J — J7l||^ + (771 — ?7l,7]2 — ??2) = 

\ ~rj2,rji —Vi + V2 = \ ||'72 — r?2||^ + (t?! - J?! i '?2 — f?2) = 

=> ||?7i — m + ?72 — ??2||^ = 
=> J71 — f?i + J72 — 572 = 

j4s can 6e uniquely decomposed into as a sum of zero, then, 

Vi ~Vi ^ 112 -V2 = 

• Let now compute 

Pho{v) 
\PhoJv)/ 

First of all, it is obvious that the constant term 770 = E(?7) and that 7712 is 
obtained by subtracting rj with all other terms of the right of the decomposition. 
Now, let us use the projector's property of embedded spaces. Indeed, as Hf C Hi, 
V i = 1, 2, it comes 

P„o{r,) = P„o{PhM) = PHfinriPQ] 

fix,) 

f is a function of Xi, so it can be decomposed into the following expression : 

<p{Xi) = v?o + ifiiiXi), ifio e Hi}, tfi £ Hi 
with ipo = E(iy3) = £(77). 

Hence, one can easily deduce Pjjo{ri), i = 1,2, as the term (pi = E(77/Xi) —£(77) 
We obtain 



(PhM\ 
Pnoiv) 
PhoM 

\PHo(r))/ 



E(77) 
E(77/Xi)-E(,7) 
E(77/X2)-E(,7) 
V77-E(7//Xi)-E(77/X2)+E(r7)/ 



(35) 



C.2. Numerical procedure 

Gauss-Seidel algorithm requires the estimation of conditional expectation at each 
iteration. To do this, we use the local polynomial estimation with a leave-one- 
out technique. To considerably reduce time cost in application, the Sherman- 
Morrison formula is exploited. 

We are going to review these methods for estimating E{Y/X = x), when Y and 
X are supposed to be real random variable. 
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The local polynomial estimation [8] consists in approximating m(x) ~ K{Y/X - 
x) by a q^^ -order polynomial fitted by a weight least squared estimation. 
An explicit solution of rh{x) is given by : 

m{x) = *(1 • • • 0)[*X(a;)i:'(a;)X(.T))]-i • 'X{x)D{x)Y (36) 

with 



(k 





D(x) 



K 



Xr, - X 



Y 



The leave-one-out technique on local estimation consists in estimating m in 
every observation point Xi, ■ ■ ■ ,Xn, i-e. computing Equation (36) when the fc*'' 
line of matrices has been removed for estimating m{Xk). It means that we would 
need to inverse *X_fc(x)I?_fc(x)X_fc(x) n times, which is very expensive. To 
avoid these expensive computations, Sherman and Morrison [24] proposed a for- 
mula : 

Lemma 2. If A is a square invertible matrix, and u, v are vectors such that 
1 + *vA-^u ^ 0, then 

+ = ^-1 _ ffl"^'' (37) 

In our problem, set Sn{x) ~ ^lL(x)D{x)lL{x). 5„(x) can be rewritten as : 

n 

Sn{x) = M^YM^) = M^YMx) +<i>k{xY<Pk{x) (38) 

i=l i^k 

V ' 

where 

$,(x) = *(i^(^V^)X,(x)), X,{x)^{l---{X,~xy), Vz = l,---,n 
h 

Thus, S-k{x), corresponding to *li{x)D{x)'K{x) when the fc*'' line has been 
removed, is of the form : 

5_fc(a;) = 5„(x) - $fc(.T)*$fc(x), V fc = I,-- - ,n 
The Sherman-Morrison formula gives : 
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As m(Xfc) = *(1 • • • 0)SzliXk) ■ *X_fc(Xfc)7^_fc(Xfc)Y, V k, it is faster to 
estimate S~^{x) and at each point of the design. 
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